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Abstract. We use a shell model to study counterion interactions in a colloidal solution. In this shell 
model, the counterions are restricted to move inside a spherical region about their host colloidal particle. 
In particular, we apply Monte Carlo simulations to derive the energy and entropy contributions of the 
effective colloidal interaction. Our result reveals an attractive electrostatic energy, which is overpowered 
by the osmotic repulsion among the counterions, as the latter can be well estimated by an ideal-gas 
approximation. We also provide an optional algorithm that enables counterion mixing between the two 
counterion clouds even when the clouds do not overlap. The residual mixing entropy of counterions gives 
a reduction in free energy that is comparable to the thermal fluctuation, suggesting a possible attractive 
mechanism between the colloidal particles under non-equilibrium condition. 

PACS. 82.70.Dd - 61.20.Ja 

1 Introduction of the DLVO theory [51[iTj]. Even after large discussions [TT1 

^^MM5S^EMS^SM, the question is still quite un- 

The phenomenon of the so-called "like-charge attraction" 

dcr debate whether the like-charge attraction origins from 

among colloidal particles has been found for almost two 

artifacts of the experimental setup. On the other hand, 

decades [1,2,3,4,5,6,7,8 in solutions of low-ionic strength. 

similar attractions have been observed in the system of 

The range of this attraction is in the order of the colloidal 

polyelectrolyte as well [2 1,2 2 , in which the phenomenon 

size, which is too large to be accounted by the van der 

is often named as "ordinary-extraordinary transition" . 

Walls interaction. Such attraction exists even when the 
counterions are monovalent, indicating a possible failure 
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Certainly counterions play a leading role in the discus- 
sions about the sources of the like-charge attraction. How- 
ever, it is important to verify if the attraction is attributed 
to their energy or entropy contribution. In a salt-free solu- 
tion, one learns from the DLVO or Debye-Hiickel[23 the- 
ory that the screening length of colloidal particles is com- 
parable with the colloidal size. This rather huge screening 
length provides a larger volume of residence for the coun- 
terions. Hence these point-like counterions about their 
host colloid can be somewhat viewed as ideal-gas parti- 
cles that are confined by the overall electric field of the 
system. Prompted by this picture, we proposed a shell 
model, in which the counterions are restricted to move in- 
side a spherical shell about their host colloidal particle. 
This shell model facilitates our task in deducing the free 
energy of two colloidal particles along with their counteri- 
ons. Our model is given in details in Sec. [3J as the results 
of simulations are presented in Sec. [4] and Sec. [5] 

The term "mixing entropy" often refers to the increase 
in entropy due to mixing of larger than two kinds of parti- 
cles. That is based on the understanding that the mixing 
of same kind of particles results in essentially no change 
of entropy in a bulk system. However, we have shown in 
our previous work[24 that there exists a non-extensive 
residual mixing entropy even if the mixing particles are 
of the same kind. This residual entropy may be negligible 
in a bulk system, but in a mesoscopic system of colloidal 
particles, the entropy effect due to the mixing of counte- 
rions results in a free energy reduce that is comparable 
with the thermal energy. In this current work, we proceed 
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with this idea of mixing, as we find from the simulation 
result in Sec. [5] that the entropy effect due to mixing of 
counterions lead to a possible like-charge attraction in the 
non-equilibrium case, where mixing is incomplete over the 
short-time scale. 

2 Attractive electrostatic energy in the DLVO 
theory 

According to the DLVO theory [MTU], the effective elec- 
trostatic interaction between a pair of colloidal particles 
is 

2 

I3F = MnR Poo ^e- K ^ 2R \ (1) 

where 

(2) 

and 

e 2 

= -r-= • 2poo (3) 
ek B T 

if the salt ions are monovalent. Also poa is the number 
density of salt ions, R is the radius of a colloidal particle, 
and d is the distance between the two colloidal particles. 
The electric potential at the colloidal surface ipo is 

according to the Debye-Hiickel theory (— Ze is the charge 
of each colloidal particle). Note that the effective energy 
described in Eq. [T] is actually a free energy (potential of 
mean force) , as it includes entropic contributions of small 
ions. To derive the pure electrostatic energy, one can apply 
the thermodynamic relation d((3F)/df3 = U. 
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For the salt-free case, is small but not exactly equal charge +e. These small counterions, modeled by point-like 
to zero due to the dissociation of water molecules. Such a particles with no size, are confined in a spherical region of 
low value of contributes to a screening length in the radius 2R about the center of the colloidal particle, as they 
order of one micrometer, as we use kR = 1 in this work, are also excluded from the colloidal volume. The reason 
Fig. [U shows the comparison of f3F, /3U and S for the case that 2R was chosen to be the size of the counterion cloud 
that Z = 100 and t = 100, as we use dimensionless units is mainly from of the fact that the size of the counterion 
such that R = 1 and fes = 1, and the reduced temperature cloud (and hence the colloidal screening length) in the salt- 
t is defined as the inverse of the Bjerrum length £b — free regime is compatible with that of a micrometer-sized 
e 2 / (AirekBT) . The entropy is derived by S — (3U — j3F, colloidal particle. We then perform Monte Carlo simula- 
with the reference point that S(d — > oo) = 0. Note that the tions, while during the course of each simulation we fix 
pure energy part /3U exhibits attractive interaction in the the distance d between the two colloidal particles and de- 
mid-range of colloidal separation, while such attractions rive the corresponding effective colloidal interaction from 
is overpowered by the entropic part, which shows osmotic the average electrostatic energy. Note that the shell model 
repulsions of the counterions. approach has been extensively studied [251126] . While the 

The DLVO expression of the effective electrostatic en- original version of the cell model allows no fluctuations 
ergy, as given in Eq.[TJ is based on the Derjaguin approximatioaf l&unterions among different cells, we use an optional 
which works only for the situation that the screening length algorithm in which the counterions are allowed to trans- 
is much smaller than the colloidal size. While this approxi- fer even if the cells are not joined together. Through this 
mation fails in the salt-free case, we perform Monte Carlo optional algorithm, along with the simplicity of the cell 
simulations to re-verify this feature of attraction in the model approach, our objective is to study the effects of 
pure-energy part in the following sections. Meanwhile, we counterion mixing in this work, 
try to examine the energy and entropy contributions with 

In this work the relevant energy is computed by con- 
more details, and check if there exists a scenario for the ef- 

sidering the electrostatic interactions only, and we use di- 

fective like-charge attractions between colloidal particles. 

mensionless variables to represent relevant physical quan- 
tities. In particular we set R= 1. Therefore the total elec- 

3 Model 

trostatic energy among the colloidal particles and their 
We set in our model that each colloidal particle has a net counterions can be written as 
charge —Ze located at its center, and the radius of a col- 
loidal particle is R. The charge of the colloidal particle is _ _ 
surrounded by Z neutralizing counterions, each bearing a ail charges y 
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where z% is the valence of the ith charge (zi = 1 for coun- tions, and for each simulation we assign a random configu- 
tcrions and zi = —Z for colloidal particles), and rij is the ration for the initial counterion distribution, in which the 
distance between the ith and jth charge. Moreover, we positions of counterions obey a uniform distribution. The 
introduce a reduced temperature t which corresponds to initialization is then followed by 10 9 simulation updates, 
the ratio of the colloidal radius over the Bjerrum length. The same procedure is repeated over 10 times to derive 
For aqueous solutions at room temperature, £-q ps 7 A, and average physical quantities. For the cases where the coun- 
therefore t is in the order of several hundreds to a thou- terion clouds do not overlap, we provide an optional mix- 
sand (since the radius of colloidal particles is in the order ing algorithm which allows counterions to jump between 
of one micrometer). The magnitude of the reduced tern- the counterion clouds, as this algorithm enables number 
perature is relatively large compared with the electrostatic fluctuations of counterions in a cloud, 
interaction between a single counterion and its host col- 
loidal particle, or the interaction between two neighboring 
counterions. For most of the time, the latter is small com- 
pared with the average thermal energy in our system with 
sparsely distributed counterions. As a consequence, the 
pairwise electrostatic interactions are rather weak, as the 
force exerted on each counterion are contributed from the 
collective charge distribution. 



4 Energetic contribution of the electrostatic 
interaction 



The result of average electrostatic energy is shown in Fig. [2] 
for the case Z = 100 and t = 100. From the squared 
data (the "non-mixing" case will be studied in a later sec- 
tion) we observe attractive electrostatic interactions as the 
In our current model we neglect the "image charge" counterion cloud starts to overlap . within the range of our 

effect mm caused by the difference in dielectric con- observation the magn i tu de of attractive energy is smaller 

stants between the colloidal particle and surrounding aque- than the thermal energy for the casc z = 1QQ ^ while 

ous solution. It has been found that such an image charge for z = 100Q thfl attractive energy becomes sign ificantly 
effect is crucial in colloidal solutions with multivalent counter^n|gj; This is in agreement with the DLV theory, al- 

Although the influence is important as well for systems though the latter ig based Qn the Derjaguin approximation 0. 

with monovalent counterions [25], we have not considered In thc next sectiong we ghall obtain the free energy of thc 

in the present work as the emphasis of the present work system and find th&t thc attractive part in e i ectrostatic 

is to bring out the importance of mixing effect. energy when the counterion clouds are over lapping will be 

To derive effective potentials and correlations, we per- overshadowed by the entropic (osmotic) repulsions, 
form Monte Carlo simulations via the Metropolis algo- The source of the attractive electrostatic energy can be 

rithm. We choose Z = 100 and Z — 1000 in our Simula- understood using the following argument: since the dis- 
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tribution of counterions is approximately spherical, the although the distribution is not identical to that of the 

net force between colloids approaches zero as the coun- DLVO result, as the decay is milder due to our sharp cut- 

terion clouds starts to separate. While the clouds over- off at the shell boundary. 

lap, the counterions inside the overlapped region are less For the case that counterion c i ouds are overlapping, 
screened and therefore can feel attractions from the col- ^he result of counterion angular distribution is represented 
loidal cores. However, this attractive energy is reduced by by a plot of counterion fraction compared with that of 
the excessive electrostatic repulsions among counterions a uniform distribution within the azimuthal angle in the 
in the overlapped region. overlapping region. From Fig. [5] one finds that when the 
To illustrate this attractive mechanism, we make a sim- count erion clouds are overlapping, some of the counterions 
pie calculation by assuming that the counterion charges diffuse out from the overlapping region, as the fraction be- 
about each host colloid are uniformly distributed. By sum- comes smaller compared with that for the non-overlapping 
ming over the Coulomb energy between all charges over cas6j the latter possessirig a uni fo rm angular distribution, 
the two colloids (except for the self energy within each The drop m number of counterions may result from os- 
colloid), one can derive motic or electrostatic repulsions of the counterions (or 

^.(2400i? 6 -6144i^+640i? 3 d 3 ^ 
90a a 

The counterion-counterion correlation function is ob- 



U HT = U{T -4 oo) = 

(6) 



where p = Z/{7 ■ 4irR 3 /3) is the homogeneous counterion tained and plotted for the non-overlapping (d = 4.1) and 

density. The behavior of U HT is shown in Fig. [3 as it re- overlapping (d = 3.0) cases in Fig. © While there seems 

veals an attractive energy of similar magnitude compared to be negligible difference between the mixing and non- 

with our simulation results mixing cases, the correlation function shows a higher peak 

Meanwhile, in our simulations we have checked about for the overlapping case compared with that for the non- 

the counterion distribution for both the non-overlapping overlapping case. From such a higher peak one may spec- 

and overlapping cases. The radial and angular distribu- ulatc that the counterions tend to 'aggregate' in the over- 

tions of counterions in the non-overlapping case are given l a PPing region, although there are chances that the result 

in Fig. H The angle represents the azimuthal angle with is due to the difference in geometric constraints between 

respect to the line joining the centers of colloidal particles, the non-overlapping and overlapping cases, 

as one finds that the counterions for the non-overlapping To verify whether there exists any attraction in the 

case are distributed uniformly angular-wise, as their dis- overall effective energy, we continue to derive the free en- 

tribution about the azimuthal angle is sinusoidal. More- ergy in the following paragraphs, as the result shows that 

over, the radial distribution exhibits a monotonic decay, the electrostatic attraction described here is indeed over- 
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powered by the entropy effect, although the result shows 
the existence of another possible mechanism of like-charge 
attraction due to counterion mixing. 
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over the electrostatic attraction, and it will be corrected 
slightly by considering the free energy change towards the 
finite-temperature regime. 



c . .. c c 5.2 Method 

5 Derivation of the free energy 

In this subsection we propose a new method to derive 

5.1 Entropy in the high-temperature limit 

the Hclmholtz free energy via Monte Carlo simulations. 
To derive the free energy one first needs the value of (3F Commonly one uses the relation 
when ft — > 0, which is equal to the negative entropy in 



d(3F 

the high-temperature limit. In our model the counteri- ®@ 



= U, (9) 

volume fixed 



ons become ideal in the high-temperature limit except for and derive the free energy by integrating the internal en- 

the restriction in their occupying volume. As the colloidal ergy over the inverse temperature /3 = 1/t. 

particles approach, the counterion clouds overlap and the We start by rewriting the partition function at inverse 

volume of occupation starts to reduce. With the use of the temperature (3o 

ideal entropy (recalling that fee = 1 in our reduced unit) 

v 

S\p =0 = 2N In V, (7) = yV^/3t/„ e -/3it/„ 

one can derive the change in entropy due to volume de- 
crease * 

where A/3 = (3q — j3\ and {) / g 1 implies the thermodynamic 

§S\p — 2N\n (^——^—^j average at the inverse temperature /3\. From Eg. [POl one 

224 \ can derive the free energy difference between /3q and /3\: 

<J}F)h-mh = -1° (tt) = -la(e-^> . (11) 
for 2 < d < 4, where d is the distance between the cen- VWi/ 

ters of the two colloidal particles (recalling that R = 1 Through the average in Eq. [11] one gains the difference 

in our model). From Fig. [7] one finds that the entropy in f3F over one simulation. This relation is exact and one 

decreases monotonically as the distance becomes smaller, does not have to apply any other integrating approxima- 

For example, 5S — —2.5836 when d = 3.5, which pro- tion. In practice, we choose (3\ = /3q/M, in which M is 

vides an increase in free energy that is crucial compared an integer (eg. M = 4 is used in our simulations). As M 

with thermal fluctuation. As we shall see next, this en- is not very large, the importance sampling still works and 

tropy decrease in the high-temperature limit is dominating one can gain the value (3F by applying Eq. [11] iteratively. 



= - 2jVln U 6 + 48d-# ) (8) 



Chi-Lun Lee: Energy and entropy effects of counterions in salt-free colloidal solutions 7 

For a system in which the dynamics may be hampered by between clouds, thus enabling counterion number fluctu- 

various local minima, this technique can be applied along ation in each cloud. Nevertheless, such diffusion takes a 

with the method of replica exchange [29 , so that free en- relatively longer time, and mixing between two counte- 

ergy is derived while frustrated dynamics can be reduced rion clouds can be rather incomplete during the observa- 

during the course of parallel simulations over various tern- tion time in a typical experiment. To estimate this effect of 

peratures. mixing, we turn our attention to the results for the mixing 

and non-mixing cases in our simulations. 

5.3 Result From Fig. [2j one finds that the electrostatic energy of 



In this work we apply the above method to derive A/3F 
between t = 100 and t = 25600, as the results for d = 3.5 
and d = 20 are shown in Table [1] Note that the resid- 
ual difference between t = 25600 and t — » oo is approxi- 
mated by (/3t/)t=25600- To compare the overall free energy, 
one has to count additionally the free energy difference in 
the high-temperature limit as derived in Sec. 15.11 For ex- 
ample, from Fig. [7] we find this entropy difference in the 
high-temperature limit to be -2.5836 between d — 3.5 and 
d = 20. Along with the result (of the mixing case) in Ta- 
ble [TJ we conclude that the overall free energy {/3F) of 
d = 3.5 is higher than that of d = 20 by the amount 
2.2805. Therefore the electrostatic attractive energy ob- 
served in Sec. 2] is heavily outweighed by the entropy ef- 
fect, the latter of which is mostly attributed to the entropy c _ i 

^mixing m 

of counterions in the high-temperature (ideal-gas) limit. 



the non-mixing case is lower than that of the mixing case. 
However, from the second law of thermodynamics [30 , one 
knows that the free energy in the mixing case must be 
lower than that in the non-mixing case, because there is 
less restriction in the mixing scenario. 

To investigate this free energy difference, we look up 
in the data in Table [1] for the case of non-overlapping 
clouds d — 20. The result shows that A(3F of the non- 
mixing case is smaller than that of the mixing case by the 
amount 0.2259. This deficit is outweighed by the entropy 
difference between the mixing and non-mixing cases in the 
high-temperature (ideal gas) limit |24] . A simple calcula- 
tion using the ideal gas entropy S = ArinV leads to a 
mixing entropy 



(2Z)l ■ 



For the case Z = 100 one obtains ^mixing = 2.8762. There- 
fore f3F of the mixing case is ultimately smaller than that 
of the non-mixing case by the amount 2.6503, which shows 
In this work we provide an optional mixing algorithm for that the free energy difference is quite compatible with the 



6 Effect of mixing 



cases that the two counterion clouds do not overlap. As to thermal energy in its magnitude. For the case Z — 1000, 
a real colloidal solution, counterions are allowed to diffuse we obtain that the free energy difference between the mix- 
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ing and non- mixing cases at d — 20 is 6.1876 using the 
same analysis. 

When the thermodynamic equilibrium is achieved, mix- 
ing of counterions is complete for all colloidal distances. 
However, the mixing between counterions of both clouds 
may be rather incomplete within the limited time of ob- 
servation at a large colloidal separation. To encourage the 
mixing process, colloidal particles will be dragged towards 
each other in the short-time interval to reduce the total 
free energy. Meanwhile, there exists a balancing mech- 
anism in which the counterion clouds tend to stay less 
overlapped to avoid osmotic repulsion. Hence our result 
gives a hint that there exists a non-equilibrium mecha- 
nism of like-charge attraction that drives the colloidal par- 
ticles towards some optimized separation, and this non- 
equilibrium mechanism is attributed to the residual mix- 
ing entropy of counterions. 

To follow the "mixing algorithm" discussed in this 
manuscript, if one considers the 'real' situation, where 
there are residual salt ions or dissociated water ions in 
the neighborhood of the counterion clouds, the counte- 
rion fluctuations of a single cloud can be easily facilitated 
through exchange with its neighboring solution, provided 
that there are ample amount of residual ions. Probably 
this explains why an effective like-charge attraction can- 
not be observed in real experiments for the case Z = 100, 
because the counterion density within the cloud is not 
much deviated with that in the ambiance. On the other 
hand, for the case Z = 1000 or even Z = 10000, the 
ionic concentration within the cloud is quite large com- 
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pared with its neighborhood. For this situation there exist 
less counterion fluctuations if the clouds do not stay close, 
since mixing cannot be easily achieved by the exchange of 
counterions of a single cloud with its surrounding, where 
residual ions are relatively sparsely distributed. 

The role of counterion fluctuations and its link to like- 
charge attractions has first been studied in the Sogami-Ise 
theory [TT], which is basically an extension of the DLVO 
theory and thus a mean-field theory, in which counterion 
fluctuations are not fully considered. Based on the DLVO 
theory they obtained an effective attraction among col- 
loidal particles via a derivation about the Gibbs free en- 
ergy instead of the Helmholtz free energy. While the coun- 
terion fluctuation considered in the Sogami-Ise theory is 
limited, due to its mean-field nature, counterion fluctua- 
tions are better described in the work of Carbajal-Tinoco 
et al. [IS] . While in the above work the authors attributed 
the like-charge attraction to charge inversions about the 
colloids due to counterion fluctuations, we have not ob- 
served in our simulations any like-charge attraction from 
the viewpoint of pure electrostatic interactions, at least in 
the current shell model description. On the other hand, 
the result of our work shows that like-charge attraction 
may arise due to entropic difference between the mixing 
and non-mixing cases, as we observe a mixing entropy even 
for identical particles (i.e, counterions). 

In the work of Polin et q/.|31j. they proposed a phe- 
nomenological expression of effective colloidal interaction 
that gives both repulsive and attractive parts. The at- 
tractive part has a longer range than the repulsive part, 
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suggesting that effective attractions emerge beyond the 
range where colloidal charges are fully screened. While the 
authors speculated that the long-range attraction is due 
to nonmonotonic distribution of counterions about their 
their host colloid, we have not derived such a distribution 
in our current simulation studies. Nevertheless, the com- 
parison between the mixing and non-mixing results also 
indicates a possible like-charge attraction whose range is 
longer than the diameter of the counterion shell, the lat- 
ter being comparable to the screening length. This is in 
agreement with the phenomenological expression by Polin 
et al.fiT\. 

7 Conclusion 

In summary, we apply our shell model to study the effects 
of counterionic energy and entropy contributions to the 
effective free energy between a pair of colloidal particles 
in the salt-free regime. We observe that the electrostatic 
energy shows an attractive part at shorter colloidal dis- 
tances, although this is overpowered by the osmotic re- 
pulsion between counterions, the latter resulting from the 
decrease of their residence volume. As our model assumes 
that counterions about each colloidal particle must reside 
in spherical shell, such a restriction may lead to a drop 
in volume that may be over estimated as the counterion 
clouds overlap. We expect the volume decrease and thus 
the entropy decrease to be smaller if the volume of counte- 
rion residence is allowed to deform in shape. Therefore we 
still cannot rule out the possibility that the electrostatic 
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attraction might eventually dominate over the osmotic re- 
pulsion. 

Meanwhile, the residual mixing entropy between the 
two counterion clouds as introduced in our previous work[Mj 
is further investigated in our current study. We argue 
that this residual entropy can lead to an effective attrac- 
tion between colloidal particles in the short-time inter- 
val, where mixing of counterions over the two colloids is 
still incomplete at larger colloidal separations. To further 
check about this idea, we suggest to study the hydrody- 
namic modes of the colloidal system using the Poisson- 
Boltzmann approach in a future work. Finally, the ar- 
gument about the residual mixing entropy is simple and 
general enough that it may be applied to explain like- 
charge attractions observed in similar systems such as 
polvelectrolvte [21H22] . 

This work was supported by the National Science Coun- 
cil of the Republic of China under Grant No. NSC-100- 
2112-M-008-005. The author would like to thank for sup- 
port from NCTS under focus group Life and Complexity. 
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tl -tl 


A/3F(d = 3.5) 


A/3F(d = 20) (mixing) 


A/3F(d = 20) (non-mixing) 


100 - 400 


-54.5479 


-54.3311 


-54.5000 


400 - 1600 


-13.4914 


-13.4268 


-13.4674 


1600 - 6400 


-3.3629 


-3.3462 


-3.3591 


6400 - 25600 


-0.8403 


-0.8366 


-0.8392 


100 - oo 


-72.5226 


-72.2195 


-72.4454 



Table 1. Free energy difference A/3F between temperatures 
ti and ti. These free energy differences are derived based on 
the algorithm as described in Sec. 15.21 For the cases where the 
counterion clouds do not overlap (i.e., d > 4) in our model, we 
provide an optional simulation algorithm in which counterions 
are allowed to jump between clouds. This jumping algorithm 
enables mixing and therefore counterion number fluctuation. 
For d — 20 we list out the results of both algorithms. 




Fig. 1. Comparison of the free energy, energy and entropy of 
the DLVO theory. In this example we use Z = 100 and t = 100, 
as the quantities are presented in dimensionless units. 
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Fig. 2. Electrostatic energy versus the colloidal distance of our 
shell model for (a) Z = 100 and (b) Z = 1000. The average 
energy is derived from Monte Carlo simulations with t = 100. 
Squares represent the results where mixing of counterions over 
the two clouds is enabled, while crosses represent the results 
of the non-mixing case. 
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Fig. 3. A mean-field approximation of the electrostatic en- 
ergy. The result is derived assuming the counterion distribution 
about each colloid is uniform. 
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Fig. 4. Counterion distributions in the non-overlapping case, 
(a) The angular distribution, (b) The radial distribution. 
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Fig. 5. The number fraction of counterions within the az- 
imuthal angle in the overlapping region for various d. The line 
corresponds to the case of a uniform angular distribution in 
each counterion cloud. 
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Fig. 6. Counterion-counterion pair distribution function g(r) 
for d = 3.0 and d = 4.1 (including both the non- mixing and 
mixing cases). 
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Fig. 7. Change of counterion entropy in the high-temperature 
limit due to volume decrease. In our shell model, the residence 
volume of counterions decreases once the two counterion clouds 
overlap, resulting in an entropy decrease. The result is derived 
using an ideal-gas law of entropy. 



